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Abstract. - We investigate the effect of diagonal disorder on bosons in an optical lattice de- 
scribed by an Anderson-Hubbard model at zero temperature. It is known that within Gutzwiller 
mean-field theory spatially resolved calculations suffer particularly from finite system sizes in the 
disordered case, while arithmetic averaging of the order parameter cannot describe the Bose glass 
phase for finite hopping J > 0. Here we present and apply a new stochastic mean- field theory 
which captures localization due to disorder, includes non-trivial dimensional effects beyond the 
mean- field scaling level and is applicable in the thermodynamic limit. In contrast to fermionic 
systems, we find the existence of a critical hopping strength, above which the system remains 
superfluid for arbitrarily strong disorder. 



P5 Ever since the seminal paper by Fisher et al. [1], the 
O disordered Bose Hubbard model has been the subject of 
, ^ t heoretical and experimental investigation. In particu- 
lar, the realization of the superfluid-Mott insulator tran- 
C^sition in a gas of ultracold bosonic atoms in an optical 
^^ lattice [2] has sparked a new wave of research on this 
^-N field. In contrast to typical solid state systems, optical 
^) lattices allow the introduction of various types of disorder 
^^in a highly controlled manner [3]. Experimentally, sev- 
T-|^eral techniques have been implemented, such as speckle 
^^ laser patterns [4-7], multichromatic lattices with non- 
OO commensurate wavelengths [8,9] or an additional species of 
^^ atoms tunneling at a considerably lower rate [10,11]. In a 
^recent experiment [6] it was possible to generate the first 
I^3D fine-grained disorder potential using a speckle laser, 
^^ providing an excellent realization of the 3D disordered 
C^ Bose-Hubbard model. Various theoretical methods have 
previously been employed to investigate the transitions be- 
tween Mott insulator (MI), Bose glass (BG) and superfluid 
(SF), such as Quantum Monte Carlo [12-16], exact diag- 
onalization [9,17-19], renormalization group [20], density 
matrix renormalization group [21] and mean- field approx- 
imations [1,22-30]. However, spatially resolved calcula- 
tions on disordered lattices suffer from finite size effects in 
the vicinity of phase borders, where the physics is dom- 
inated by rare events, while an arithmetically averaged 




Fig. 1: Within SMFT, the multiple site lattice model is ap- 
proximated by an effective single-site problem, where a site is 
coupled to a bath of mean- field parameters (MFPs). Disorder- 
induced fluctuations of the MFPs are accounted for by a sta- 
tistical distribution P{ip). 

mean-field theory is incapable of describing the Bose glass 
phase at any finite hopping amplitude and T = [29]. 

Ultracold bosonic atoms in a sufficiently deep optical 
lattice at moderate filling are well described by the single 
band Bose-Hubbard (BH) Hamiltonian [31] 
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where h\ is the bosonic creation operator for an atom in 
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the lowest Wannier state at lattice site z, J is the hop- 
ping amplitude, U describes a short-ranged interaction, 
(i, j) denotes the sum over all pairs of neighboring sites 
and /i is the chemical potential. In performing the low- 
est band approximation, it is implicitly assumed that the 
temperature is sufficiently low to suppress all contribu- 
tions of states in higher bands. Since all system prop- 
erties only depend on the ratios of the energy scales, we 
choose to work in units of U = 1. Diagonal disorder is 
parametrized by on-site energies e^, which we will assume 
to be independently and identically distributed ^ accord- 
ing to a box distribution p{e) = 6(A/2 — |e|)/A with the 
disorder strength parameter A, while the hopping ampli- 
tude J is assumed to be site-independent. On a mean- field 
level, the superfluid-insulator transition can be captured 
by the variational bosonic Gutzwiller ansatz [22] 
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which is a good approximation in high dimensions (coor- 
dination number Z) and delivers the exact ground state 
in the weak tunneling limit J ^ 0, as well as in the non- 
interacting limit U ^ [32] . The mean-field ground state 
is determined by minimizing the energy expectation value 
(GW|7Ybh|GW) with respect to all the amplitudes {fn } 
under the constraint (GW|GW) = 1. For a pure system 
(A = 0) at T = this leads to the same ground state as 
that of the site decoupled mean- field Hamiltonian [27] 
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where the mean- field parameters (MFPs) '^i = (bi) are de- 
termined self-consistently. In the ground state, all MFPs 
have the same complex phase and, due to the global U{1)- 
symmetry of TYbh, can be chosen to be real and non- 
negative. 

For A = a vanishing order parameter (defined as 
any of the identical MFPs) indicates an insulating state 
of the system, while a finite value indicates a Bose con- 
densed (SF) phase. The spatially resolved, self-consistent 
Gutzwiller approach has been used to study disordered 
bosons in optical lattices within finite size simulations 
[24,25], but suffers of the disadvantage that it overesti- 
mates phase coherence, especially in the SF phase in the 
vicinity of the phase borders and cannot describe the BG 
phase at T = 0. However, simulating disorder effects in 
finite systems is a delicate problem, since rare events may 
strongly influence physical observables such as the excita- 
tion spectrum. 

Here we present and apply a stochastic mean-field the- 
ory (SMFT) for disordered bosons, which extends the 




■"^This is a justifiable approximation considering new experimental 
techniques [6] 



Fig. 2: Results from a spatially resolved, bosonic Gutzwiller 
calculation (27 000 sites) in the deep SF regime at J = 0.067, 
Z = 6,/i=l, A = 0.5, where this method becomes exact 
(Gross-Pitaevskii regime with depletion). The left figure is 
a color-coded histogram (corresponding to the 2D probabil- 
ity density function), showing the strong correlation between 
the on-site energy e and the on-site MFP ip (correlation coeffi- 
cient -0.9964), which is not neglected within SMFT. The cen- 
tral figure is a color-coded histogram of the on-site energy and 
nearest neighbor MFPs, which are only weakly correlated (cor- 
relation coefficient -0.0349) and neglected within SMFT. The 
right figure is a comparison of a self-consistently determined 
distribution from SMFT for these parameters, compared to a 
ID histogram of MFPs from the spatially resolved Gutzwiller 
calculation (fluctuations originate from the finite system size). 
Excellent agreement is achieved in this limit. 



self-consistent Gutzwiller approach to the thermodynamic 
limit and is free of finite-size effects. Its numerical effi- 
ciency and absence of finite size effects make this method 
a good candidate for future analysis of more complex sys- 
tems, such as multicomponent gases in disordered lattices. 
This is achieved by a probabilistic description of an infi- 
nite system, by using a probability density function (PDF) 
P(V^) to allow disorder-induced fluctuations of the MFPs. 
Probabilistic descriptions have been successfully applied 
to disordered systems previously, such as by [33] for anti- 
ferromagnets. 

The many-particle, multiple-site Bose-Hubbard model 
is thereby reduced to an effective single-site mean- 
field problem, the solution of which entails the self- 
consistent determination of P('0). To derive the SMFT 
self-consistency condition for P('0), we pursue the self- 
consistent mean-field approach and assume that the on- 
site energy of a given site and the MFPs of the neigh- 
boring sites are uncorrelated. To justify this approxima- 
tion, we performed a spatially resolved Gutzwiller calcula- 
tion, which is known to give a good approximation for the 
ground state for weak interactions U <^ JZ. In Fig. 2 we 
plot a histogram, showing that the correlation between an 
on-site energy and the nearest neighbor MFP is weak and 
the above approximation is justified for weak interactions. 

Considering an arbitrary site i with energy e^, the 
further quantity determining its mean-field ground state 
is the scaled sum of MFPs from the neighboring sites 
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T] = J ^^ n 7=1 ^j distributed according to the PDF 
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Hence, once the self-consistent solution P('0) is known, 
any disorder averaged expectation value of a single-site 
operator can be expressed as 

T^) = j dep{e) J dn Qiv) (gs(e, r,) |i|gs(e, rj)), (5) 

where |gs(e, 7^)) is the ground state of 

n = r]{h^ ^h)^{e-ii)h^h^-h^h^hh. (6) 

The self-consistency condition requires that if the on- 
site energy e is randomly chosen from p(e) and Z MFPs 
are drawn from P('0) to account for an effective tunneling 
from the nearest neighbors (or equivalently 77 is drawn 
from Q (/?)), the calculated expectation values {h) have 
to be distributed according to the initially assumed PDF 
P('0) (illustrated in Fig.l). This can be expressed by the 
self-consistency equation 
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is the conditional probability density for a site having the 
MFP tp if the external coupling rj is given and the disorder 
energy is distributed according to p{e). 

As a next step it is of interest to characterize the dif- 
ferent phases once the distribution P('0) has been deter- 
mined. The condensate fraction within the SMFT is given 

by 

/c = {b)ym, (8) 

which reveals that the system is in the superfluid phase as 
soon as P{tp) 7^ S{iIj) (i.e. /c > 0). Rigorously speaking, 
this should be referred to as the condensed phase, as a 
finite value of the average MFP indicates the existence of 
a macroscopically occupied single-particle state. To dis- 
tinguish the MI from the BG phase, a further quantity, 

such as the compressibility k, = -^ (which vanishes only 
in the MI), has to be considered. 

For general parameters /i, J, Z, and A, the distribu- 
tion P('0) fulfilling the self-consistency condition (7) can- 
not be determined analytically. We determine P(V^) by a 
numerical iterative procedure, beginning with any distri- 
bution, other than the insulating solution Pmi('^) = ^(V^), 
which always fulfills the self-consistency equation as it is 
a fixed point of the iterative mapping. We have verified 
numerically that there always exists a unique attractive 
self-consistent solution, which we identify as the physical 
distribution (which also minimizes the grand canonical po- 
tential). A multiple grid discretization procedure, yielding 




Fig. 3: Self-consistently determined MFP probability density 
distributions at fixed /i= 1.0, J = 0.05, Z = 6 and increas- 
ing disorder strength A. The distributions are normalized to 
their maximum value for visual clarity. For A ^ the PDF 
converges to a shifted ^-distribution, recovering to the well- 
known usual MFT for the pure system. With increasing A the 
disorder induces fluctuations in the MFPs, i.e. P{ip) broadens. 



a high resolution to detect the superfluid-insulator transi- 
tion at small values of ?/^ was used for the numerical tab- 
ulation of P('0). 

Typical results for such distributions are shown in Fig. 3 
for a variety of increasing disorder strength values A. At 
A = the distribution P(V^) is a ^-function at the value 
of tjj corresponding to the solution of the usual Gutzwiller 
bosonic MFT (in the disorder- free case). In the presence 
of disorder (A > 0), P('0) acquires a finite width in the 
SF phase. This can be understood to have two origins: 
Fluctuations of the on-site energy e necessarily lead to a 
variation in the calculated MFP (b) (for non-zero MFPs 
from the neighboring sites). Subsequent additional fluctu- 
ations in the MFPs furthermore enhance the fluctuations 
of (b). By decreasing the hopping strength for fixed dis- 
order, we find that the system is always driven into an 
insulating state with ip = 0. Care has to be taken at 
the SF/BG transition, where it has to be ensured that the 
distribution is independent of the numerical discretization. 
Since no averaging of the MFPs is performed, the results 
of the SMFT depend nontrivially (beyond JZ-scaling) on 
the dimensionality of the system: with increasing dimen- 
sionality the BG region in the phase diagram gives way to 
the superfluid phase. Using the properties of the convolu- 
tion, it can be seen that in the limit of infinite dimensions 
the arithmetically averaged MFT is recovered, where the 
BG can only exist at J = 0. In this sense, the existence 
of the BG depends on fluctuations of the MFPs, which 
are accounted for by the probabilistic description within 
the SMFT. Since the relative fluctuations of the MFPs 
decrease with the dimension of the system, the BG phase 
exists in a larger region of the phase diagram (scaled with 
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disorder averaged mean field parameter \|/ 
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Fig. 4: Plots of the arithmetic mean order parameter ip — f dip P{'iIj) ip (upper row) and the compressibility hi (lower row) show 
the stochastic mean- field phase diagram for three different disorder strengths and Z = 6. The white lines indicate the phase 
borders and black corresponds to the value zero for all plots. For weak disorder (A = 0.1) the phase diagram closely resembles 
that of a pure system. With increasing disorder the Mott insulating regions (V^ = 0, a^: = 0) shrink and are completely replaced 
by the Bose glass (-0 = 0, Ai: > 0) and superfluid ('0 > 0, At > 0) phases at A = 1. The areas for /i < with '0 = and n — ^ 
correspond to the vacuum, (color online) 



JZ) in low dimensions. Furthermore, the question of a 
critical dimension for the crossover from a direct transition 
between the MI and SF and an indirect transition, always 
occurring via the BG phase, arises^. The latter scenario 
is well-established in one and two dimensions, while in the 
limit of high dimensions, where SMFT is known to become 
exact, a direct transition at the tip of the Mott lobes is 
predicted. 

The resulting phase diagrams for a 3D lattice are shown 
in Fig. 4. The transitions between the MI, BG and SF 
phases are of second order, since the density (n(/i)) varies 
continuously, but the compressibility is discontinuous at 
the transition points, as shown in Fig. 5. In the BG phase, 
the compressibility is proportional to the density of sites at 
positive integer values of the effective chemical potential 
ji' = ji — e. Therefore it takes on constant positive values 
on finite intervals of /i in the BG phase for a box disorder 
distribution, while it varies continuously within the SF 
phase, as can be seen in Fig. 5. 

The SMFT predicts that for every fixed value of /i, there 
exists a certain value JZc{/J^) above which the system is 
always in a superfluid state, independent of the disorder 
strength A. Typical phase diagrams in the JZ-A-plane for 
constant /i = 0.4 and for integer filling n = 1 are shown 
in Fig. 6. With increasing disorder strength at constant /i. 



- - JZ = 0.009 — JZ = 0.099 — JZ = 0.219 




^After completing this work, a preprint also addressing this ques- 
tion appeared [34] 



Fig. 5: The compressibility hi — ^ for intermediate disorder 
strength A = 0.6 and Z = 6 as a function of /x. For small 
hopping parameters (blue line JZ = 0.009) the system is driven 
through an alternating sequence of MI and BG phases. For 
stronger hopping (red line JZ = 0.099) the system enters and 
remains in the superfluid phase from some value of /i onwards, 
where it varies smoothly (color online). 



we find that the superfluid-insulating phase border moves 
to smaller values of JZ, while fluctuating in A with a 
periodicity of 2U . This periodic behavior at small fi/U 
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and upcoming experiments on disordered bosons in optical 
lattices. 

Being a single-site theory, the SMFT is numerically 
less demanding than other methods describing the BG 
transition, such as Quantum Monte Carlo calculations or 
renormalization group studies. This allows an extension 
of the method to more complex systems, such as multi- 
component mixtures. However, the single-site nature has 
the drawback that spatial information is lost. The stochas- 
tic approach chosen here may furthermore be extended in 
future to describe fluctuations of different origin, such as 
thermal or quantum fluctuations. 



Fig. 6: SMFT phase diagram in the JZ — A-plane for fixed 
/x = 0.4 (gray line, small black dots) and fixed density n — \ 
(blue line and larger circles, with error ±0.003 in JZ\ Inset: 
illustration of the bosonic background screening the strongly 
disordered spatial potential. The MI/BG phase borders occur 
at constant A, since the mean- field state cannot depend on JZ 
if all MFPs -0 vanish, (color online) 



originates from the relative weight of Mott insulating lobes 
within the disorder interval. 

This can be understood qualitatively by interpreting the 
disordered BH model in terms of the pure phase diagram, 
where the disorder corresponds to a whole range of dif- 
ferent values of /i contributing to each point in the disor- 
dered phase diagram. Deep wells (i.e. lattice sites with 
low effective chemical potential) are successively filled up 
with a suitable number of particles, forming a bosonic sea 
which effectively screens the disordered potential - on top 
of which it is energetically favorable for the remaining par- 
ticles to delocalize (illustrated in the inset in Fig. 6). In 
the limit of strong disorder, it should be kept in mind 
that these predictions are only valid in a regime where the 
single band Bose-Hubbard model is justified. 

In conclusion, we have developed a stochastic mean-field 
approach to the disordered Bose-Hubbard model in three 
spatial dimensions. By working with the full distribution 
function of MFPs without averaging these, we are able to 
describe the Bose glass phase and the underlying localiza- 
tion of bosons within a local approach at T = 0, which 
becomes rigorous in the limit of high spatial dimensions. 
In contrast to spatially resolved Gutzwiller calculations, 
which predict condensation into a superposition of dis- 
tant localized single-particle states (i.e. overestimate the 
coherence) in the BG regime, the SMFT does not suffer 
from this problem. Furthermore, we observe a direct tran- 
sition between Mott insulator and superfluid in the pres- 
ence of disorder and find that superfluidity persists above 
a critical hopping strength for arbitrarily strong disorder 
at fixed /i, due to screening of strong potential fluctuations 
by accumulation of bosons. These flndings and quantita- 
tive predictions are of immediate relevance for current [6] 
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